Kinetics of Phase Separation in Fluids: A Molecular Dynamics Study 
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We present results Irom extensive 3-d molecular dynamics (MD) simulations of phase separation 
kinetics in fluids. A coarse-graining procedure is used to obtain state-of-the-art MD results. We 
observe an extended period of temporally linear growth in the viscous hydrodynamic regime. The 
morphological similarity of coarsening in fluids and solids is also quantified. The velocity fleld is 
characterized by the presence of monopole-like defects, which yield a generalized Porod tail in the 
corresponding structure factor. 
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The nonequilibrium evolution of a phase-separating 
binary mixture, A-l-B, is a coniplex nonlinear process 
[l|. This problem has attracted much research inter- 
est both computationally [2] and experimentally [3|. 
The growth of A-rich and B-rich domains during phase 
separation is a scaling phenomenon. The two-point 
equal-time correlation function, C^^{r,t), which char- 
acterizes the domain morphology and growth, scales as 
C^^{r,t) — g{r /(.{€)) [J|. Here, g{x) is a scaling func- 
tion independent of time. The average domain size (.{t) 
grows with time t as 



i{t)^e. 



(1) 



The growth exponent a depends upon the transport 
mechanism which drives segregation. For diffusive dy- 
namics, i ^ t^'^, which is referred to as the Lifshitz- 
Slyozov (LS) law [IJ. The LS behavior is the only 
growth law expected for phase-separating solid mix- 
tures. However, for fluids and polymers, one expects 
faster growth at large length scales where hydrodynamic 
effects are dominant. For d — 5, convective transport 
yields additional growth regimes [5] with 



a = 1, i{t) < 4 
a = 2/3, e{t) » 4 



(2) 



In Eq. ([2]), the inertial length £i^ [~ rf /{frf), rj, p and 7 
being the shear viscosity, density and interfacial tension] 
marks the crossover from a low-Reynolds-number vis- 
cous hydrodynamic regime to an inertial regime. There 
has been experimental evidence [6| for a crossover from 
diffusive to viscous growth. However, no experimental 
observation of an inertial regime has been reported. 

While recent focus has turned to systems with real- 
istic interactions and boundary conditions [2, y], our 
understanding of segregation kinetics in bulk fluids re- 
mains far from complete. The viscous regime has been 
observed in numerical studies using the phenomenolog- 
ical Model H [TJ, |8[. Further, both viscous and inertial 



regimes have been observed in lattice Boltzmann sim- 
ulations [Sl, [lOJ. However, molecular dynamics (MD) 
methods, where hydrodynamics is automatically in- 
built, have rarely been used to study domain growth, 
primarily due to heavy computational requirements. To 
the best of our knowledge, the first MD study was by 
Ma et al. [ll| , who did not find a signature of viscous 



growth. In a later MD simulation, Laradji et al. [12 1 
observed linear domain growth over a small interval in a 
binary Lennard- Jones (LJ) fluid. More recently, Thakre 
et al. [13l| used a similar model to study the crossover 
from diffusive to viscous dynamics. However, they do 
not observe linear growth in the post-crossover regime. 
In all these cases, MD results have been obtained for 
low-density fluids over very limited time-windows, and 
conclusions drawn from these should not be taken se- 



riously. In related work, Kabrede and Hentschke [14 1 
found a ~ 0.5 in MD simulations of gas-liquid phase 
separation. In this letter, we present results from large- 
scale MD simulations in conjunction with a numeri- 
cal renormalization-group (RG) procedure [l5|. These 
state-of-the-art results will serve as a valuable reference 
for experimentalists and theorists in this area. 

Following Das et al. [l6|, we have employed a sym- 
metric model where particles of diameter a interact 
via the potential V{rij) = U{rij) — U{rc) — {rij — 
rc)dU{rij)/dr^j\r^.^r,- Here U{rij) = 4eap[{cr/rij)^'^ - 
{(y/rij)^] is the LJ potential; and rij = |ri — rj\, Vc = 
2.5a and a, /3 = A, B. We chose eaa = ebb = 2eAB — e, 
so that phase separation is favored energetically. All 
particles were assigned equal mass m. We set m, a, e 
and ks to unity. An incompressible fluid (p — 1) is 
studied for which phase separation sets in at a critical 
temperature Tc ~ 1.423 [l6j, well separated from the 
gas-liquid and liquid-solid transitions. A total number 
of 262144 particles were confined in a cubic box of size 
64'^ with periodic boundary conditions in all directions. 
The MD runs were performed using the standard Ver- 
let velocity algorithm [l7[ with a time step At — O.OIt, 
where the LJ time unit r = (mu^/e)^" = 1. The tem- 



perature T was controlled by a Nose-Hoover thermostat 
(NHT) 17], which is known to preserve hydrodynam- 
ics. Homogeneous initial configurations were prepared 
by equilibrating the system at T = 10. At t = 0, the 
system is quenched to T < Tc- All the results pre- 
sented here were obtained by averaging over 5 indepen- 
dent runs at a quench temperature T — O.llTc- 



t = 




FIG. 1: Evolution snapshots for a 50:50 binary Lennard- 
Jones fluid after quenching from the high-temperature ho- 
mogeneous phase to T = 0.77Tc. The A and B particles are 
marked in black and grey, respectively. 

In Fig. [1] we present evolution pictures at different 
times. As expected for a symmetric (critical) compo- 
sition, a bicontinuous domain structure is seen. The 
snapshot at t = corresponds to the homogeneous 
state immediately after the quench. The snapshot at 
t = 7000 corresponds to the situation where the sys- 
tem has almost completely phase-separated. Note that 
this time interval is more than an order of magnitude 
larger than earlier MD studies. While domains grow 
without encountering any perceptible size effects within 
this time-window, finite-size effects are seen beyond it 

To characterize the domain morphology, we calculate 
the correlation function as 



C^^{r,t) 






(3) 



where the order parameter ^(r, t) [= XAif, t) — XBir, t)] 
is the local concentration difference between A and B 
species. The angular brackets in Ec^. ([3]) denote sta- 
tistical averaging. We use a coarse-graining (numerical 
RG) procedure [l5[ to obtain the pure domain structure 
by eliminating thermal fiuctuations in the snapshots of 
Fig. [TJ Fig. [5] shows the scaling plot of C^^{r,t) vs. 



r/l. The average domain size £ is defined as the first 
zero crossing of C^^ (r, t) , which is computed from the 
coarse-grained order parameter. Our correlation func- 
tion data is comparable in quality to that obtained 
from a Langevin simulation IS|. The neat data collapse 
over an extended interval shows that a scaling regime 
is reached. 




FIG. 2: Scaling plot of the correlation function, C^^{r,t) 
vs. r/£, for 3 different times. The inset shows the scaling 
plot of the structure factor S^^{k,t) for the same times. 
The solid and dashed lines denote analogous data from the 
Kawasaki-Ising model. 

We would like to make a quantitative comparison of 
our MD results for C^^ (r, t) with those for segregation 
in the Kawasaki-Ising model (KIM), where the struc- 
ture and dynamics are much better understood. The 
lines in Fig. [2] denote C^^{r,t) vs. r/i, obtained from 
a Monte Carlo (MC) simulation of the KIM ^ with 
critical (50:50) composition. The MC scaling function 
is in excellent agreement with our MD data. In the 
inset of Fig. [51 we show the corresponding results for 
the scaled structure factor, Sjf,^{k,t)£~^ vs. k£, where 
S^^{k,t) is the Fourier transform of C^^{r,t). Again, 
a good data collapse is obtained confirming the scaling 
form, S,p,p{k,t) — i'^f(ki) [ij. The agreement with the 
corresponding KIM result is demonstrated again in the 
inset. Here, the decay of the tail with a power law, 
SjMj ~ k~'^, is consistent with the expected Porod's law 
[24 [22!, S{k,t) ~ yfc-(''+"), for ordering dynamics in 
d = 3 with a scalar order parameter (n = 1). (In an ex- 
tended publication, we will present results for morpho- 
logical characteristics like the Tomita sum rule, Yeung- 
Furukawa law for fc — > 0, etc.) The excellent agree- 
ment of the scaling functions confirms a close similarity 
of structures formed during phase separation in fiuids 
with those for solid mixtures. In a related context, Puri 
et al. 23^ have emphasized that domain growth mor- 



phologies are approximately independent of the kinetic 
mechanism of coarsening. This should be contrasted 
with a Cell Dynamical Systems study by Shinozaki and 
Oono [8|, who argued that there were different scaling 
functions for phase separation in alloys and fluids. 
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FIG. 3: Plot of average domain size, lit'), as a function of 
time t. The insets show the instantaneous exponent ai vs. 
1/i? (upper) and ai vs. 1// (lower). A detailed explanation 
of the insets is provided in the text. 

Next, we focus on the time-dependence of the domain 
size. In Fig. [3l we plot i vs. t. The growth at later 
times (t > 2000) is clearly linear, but the earlier-time 
data deviates somewhat. In fact, a least-squares fit to 
the form i(t) = B + At"" in the range t e [0, 2000] gives 
an exponent a ~ 0.75. Essentially, the diffusive regime 
(with £ ^ ti^/"^) is very short-lived at this temperature. 
A gradual crossover to the viscous regime (with £ ~ 
t) starts very early, extending over a large fraction of 
the time window, where the effective exponent is rather 
high. 

The top curve in the upper inset of Fig. |3] shows 
the instantaneous exponent ai[— d(ln£)/d(lnt)] vs. l/£. 
Because of the significantly large off-set value in £ as 
i — ^ 0, ai tends to its asymptotic value in a linear fash- 
ion only in the limit ^ — > oo [18]. Without a knowl- 
edge of the functional form of Ui for the whole range 
of l/£, the dashed straight line serves as a guide to the 
eye and suggests at -^ 1. The lower curve in this in- 
set is obtained from the coarsening dynamics when an 
Andersen thermostat (AT) is used instead of an NHT. 
The AT, where T is controlled by letting the particles 
collide randomly with a heat bath, is stochastic in na- 
ture and does not model hydrodynamics. Thus, it is 
expected to provide a diffusive growth of domains with 
a = 1/3. Indeed our result is consistent with that ex- 



pectation, as shown by the dashed straight line which 
extrapolates to a = 1/3 in the upper inset. This con- 
firms the utility of the NHT in studying hydrodynamic 
phenomena in domain growth. We also emphasize that 
MD with an AT offers a more realistic way of modelling 
diffusive phase separation than the commonly-used but 
unphysical KIM. 

There is another instructive way of investigating at. 
In hindsight, we introduce a time tg and assume that 
segregation kinetics follows a power-law behavior with 
time t' — t — to'. 



£'{t')^£{t)-£{to)^At" 



(4) 



Then, we calculate the exponent ai — d{ln£')/d{\nt'). 
For linear growth, Eq. (j3]) is invariant under an arbi- 
trary choice of tp. Thus, if to is chosen appropriately, 
tti ~ 1 for all values of t' . However, as noted by other 
authors [8|,ll8|, in computer simulations of finite systems 
one finds an oscillation of a^ as £ ^- cxd, with growing 
amplitude around the expected value. This is due to 
increasing separation between the domains of like par- 
ticles, thus delaying collisions between domains of large 
size. The lower inset of Fig. [3] plots ai vs. l/£' for 
to — 2500, which lies in the linear region. Indeed, this 
plot is consistent with the above expectation, and ai 
oscillates around the mean value a — 1. 
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FIG. 4: Pattern formation in the velocity field a.t t = 7000. 
(a) Projection of velocity field onto a 2-d cross-section, (b) 
Same as (a) but for the coarse-grained velocity field, (c) 
Plot of correlation function, Cvv{r,t) vs. r, for the velocity 
field, (d) Plot of corresponding structure factor, Svv{k,t) 
vs. k, on a log-log scale. 

Finally, we turn our attention to the pattern dy- 
namics of the velocity field. Fig. SJ^a) shows a 2-d 
cross-section of the system at t = 7000, with particle 
velocities being projected onto this plane. While the 
orientations of velocity vectors look fairly random at 
the microscopic scale, structure starts emerging upon 



coarse-graining over larger length scales, as seen in Fig. 
lUJb). In Figs. \M,c) andSJcl), we present plots of the 
velocity correlation function Cw [f, t) and structure fac- 
tor Svv{k,t), analogous to the density field. The decay 
of the structure factor tail as Sw ~ k~^ is consistent 
with the generalized Porod's law 22] for ordering of a 
3-component (n — 3) vector field in d = 3. This is 
indicative of pattern formation with monopolc-like de- 
fects, which can be seen in the coarse-grained snapshots 
of the velocity field. To study growth in this ordering, 
as seen in the lattice Boltzmann simulations [9|, |lO[ , we 
need coarse-graining over a larger length scale which is 
not accessible with the present system size. 

In summary, our results from extensive MD simula- 
tions of a binary L J fluid unambiguously confirm the lin- 
ear growth law in the viscous hydrodynamic regime. We 
use a numerical RG technique to obtain noise-free data 
from our MD studies. Even though the growth mech- 
anisms in fluids are different from solids, the domain 
morphologies are comparable in the two cases. Similar 
studies at higher quench temperatures will be useful to 
identify diffusion-driven growth and crossovers at early 
times. Note that while our choice of system size was 
appropriate for understanding viscous growth, this size 
is not large enough to study dynamics in the inertial 
regime. To achieve the latter goal, sophisticated mod- 
elling at the multi-scale level is required, in addition to 
usage of parallel programming and a graphics card. 
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